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ABSTRACT 

We address the highly debated issue of constraining the 7 -ray emission region in blazars 
from cross-correlation analysis using discrete correlation function between radio and 7 -ray 
light curves. The significance of the correlations is evaluated using two different approaches: 
simulating light curves and mixed source correlations. The cross-correlation analysis yielded 
26 sources with significant correlations. In most of the sources, the 7 -ray peaks lead the radio 
with time lags in the range h- 20 and h- 690 days, whereas in sources 1633H-382 and 3C 345 
we find the radio emission to lead the 7 rays by -15 and -40 days, respectively. Apart from 
the individual source study, we stacked the correlations of all sources and also those based 
on sub-samples. The time lag from the stacked correlation is h- 80 days for the whole sample 
and the distance travelled by the emission region corresponds to 7 pc. We also compared 
the start times of activity in radio and 7 rays of the correlated flares using Bayesian block 
representation. This shows that most of the flares at both wavebands start at almost the same 
time, implying a co-spatial origin of the activity. The correlated sources show more flares and 
are brighter in both bands than the uncorrelated ones. 

Key words: galaxies: active - galaxies: jets - galaxies: nuclei - gamma rays: galaxies - radio 
continuum: galaxies - radiation mechanisms: non-thermal 


1 INTRODUCTION 

Blazars (i.e., flat-spectrum radio quasars, FSRQs, and BL Lac ob¬ 
jects) are a type of active galactic nuclei (AGN) with a jet oriented 
close to the line of sight. The broadband spectral energy distribution 
of blazars is characterized by the presence of two broad bumps. The 
low energy component (from radio to UV or, in some cases. X-rays) 
is produced via synchrotron radiation by relativistic electrons in the 
jet plasma. The high energy component, generally peaking at the 7 - 
ray regime, can be produced by inverse Compton (IC) scattering of 
the same electrons responsible for the synchrotron emission (lep- 
tonic models) or emission resulting from the cascades initiated by 
photo-pair and photo-pion production of ultrarelativistic hadrons in 
the jet (hadronic models). For further information on leptonic and 
hadronic blazar models, see Bottcher et al. (2013). 

Despite many studies, the knowledge on the location of the 
7 -ray emission and its mechanism remains uncertain. Several mod¬ 
els have been proposed regarding the 7 -ray emission site relative to 
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the central engine in blazars. Some of them constrain the location 
closer to the supermassive black hole (< 0 . 1-1 pc) based on the ob¬ 
served rapid 7 -ray variability in few sources (e.g., Tavecchio et al. 
2010; Foschini et al. 2011). In such cases, the 7 -ray emission 
could be generated by IC process on the external photons (EC), 
in which the soft photons are directly from the accretion disc (e.g., 
Dermer & Schlickeiser 2002) or from the broad line region (BLR; 
e.g., Sikora et al. 1994). 

Multifrequency and relative timing analysis of individual 
sources constrains the 7 -ray emission site to be on parsec scales 
in the relativistic jets of blazars. Results from single-dish ra¬ 
dio/mm observations shows that the 7 -ray flares are typically 
found between the onset and the peak of the mm outburts 
(Lahteenmaki & Valtaoja 2003; Leon-Tavares et al. 2011). Con¬ 
nection of 7 -ray flares to the propagation of superluminal compo¬ 
nents in the jets observed using very long baseline interferometry 
(VLBI) suggests that the kinematics of shocks or shock-shock in¬ 
teraction at parsecs down the jet as possible sites for the 7 -ray emis¬ 
sion (e.g., Jorstad et al. 2001, 2013; Agudo et al. 201 la, b). Under 
this scenario, the 7 -ray emission could either be generated via IC 
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scattering of seed photons that originate as synchrotron emission in 
the jet (synchrotron self-Compton or SSC; e.g., Bloom & Marscher 
1996) or from the jet consisting of highly-relativistic spine sur¬ 
rounded by a mildly-relativistic sheath as a source of EC (spine- 
sheath model; Marscher et al. 2010). 

Thus, one way to address this question is by multifrequency 
correlation studies that can constrain the location of emission re¬ 
gion. Recent works such as, Cohen et al. (2014), Fuhrmann et al. 
(2014) and Max-Moerbeck et al. (2014a), investigate the correla¬ 
tions between 7 rays and optical and radio light curves. These 
works have utilised the 7 -ray data provided by the ongoing 
successful mission of Fermi Gamma-ray Space Telescope along 
with its primary scientific instrument the Large Area Telescope 
(Terrai/LAT) providing unprecedented coverage of the 7 -ray sky. 

In this work, we discuss the results from the correlation anal¬ 
ysis between the 37 GHz radio light curves obtained from the 
Metsiihovi AGN monitoring programme and the Fermi 7 -ray light 
curves of 55 blazars. The paper is organised as follows. In Section 2 
we describe the sample definition used for selecting 55 blazars. We 
describe the observation and data reduction methods in Section 3. 
Sections 4 and 5 describe the correlation and various significance 
methods used in this work. The results and discussion of the cor¬ 
relations are presented in Sections 6 and 7 before summarising the 
work in Section 8 . 

We assume a A cold dark matter cosmology throughout this 
work with Flo = 68 km s“^Mpc“^, Gm = 0.3, and Ga = 0.7 
(Planck Collaboration et al. 2014). The term radio core used in 
many places in this work refers to the compact unresolved feature 
in the VLBI maps that is usually located at one end of the jet. 


2 SAMPLE DEFINITION 

For this work, we considered all AGNs from the Metsahovi mon¬ 
itoring programme with observations since 2008.6 that are also 
part of the second Fermi Gamma-ray catalogue (hereafter 2FGL; 
Nolan et al. 2012). We required the 37 GHz mean flux density to 
be above 1 Jy for the time period 2008.6-2013.6. We obtained 55 
AGNs which were also variable at the 3cr level according to a 
test. Our sample could be sub-divided into 40 FSRQs, 14 BL Lacs 
and 1 radio galaxy. The source names, optical classification and 
redshifts are given in Table 1 . 

3 OBSERVATIONS AND DATA REDUCTION 
3.1 FermifLKT 

The FermifFPS is an electron-positron pair conversion tele¬ 
scope covering the energy range from 20 MeV to >300 GeV 
(Atwood et al. 2009). It operates in a survey mode observing the 
entire sky every 3 hr. This makes Fermi an ideal source of 7 -ray 
data for our analysis. 

The 7 -ray fluxes for the energy range 0.1-200 GeV were 
obtained by analysing the FermilLKI data from 2008 August 4 
to 2013 August 3 (encompassing 5 yr of 7 -ray data) using the 
Fermi Science Tools ^ version v9r32p5. Following the data selec¬ 
tion recommendations^, we select photons in event class 2 , with 

^ http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone 
^ http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/ 
Cicerone_Data_Exploration/Data_preparation.html 


a further selection on the zenith angle > 100 ° to avoid con¬ 
tamination from the photons coming from the Earth’s limb. The 
photons were extracted from a circular region centered on the 
source, within a radius of 15°. The instrument response functions 
P7REP_SOURCE_V15 were used. 

We implemented an unbinned likelihood analysis using gtlike 
(Cash 1979; Mattox et al. 1996), modelling all the sources within 
15° from our source position (region-of-interest, ROI) obtained 
from the 2FGL. We fixed the model parameters of sources with 
test statistic (TS) < 4 to the 2FGL value and also for sources 
> 10° from ROI. The normalization was left free, while the spec¬ 
tral indices for all the sources were fixed at the 2FGL value. The 
Galactic diffuse emission and the isotropic background (sum of 
extragalactic diffuse and residual instrumental backgrounds) were 
also modelled at this stage using the templates ’gll_iem_v05.fits’ 
and ’iso_source_v05.txt’ provided with the Science Tools. Our final 
fluxes were obtained from 7 d integrations with a detection criterion 
such that, the maximum-likelihood TS (Mattox et al. 1996) exceeds 
four (~ 2a). Bins with TS < 4 or when the predicted number of 
photons are less than four, 2a upper limits were estimated using the 
profile likelihood method (Rolke et al. 2005). 

In addition to the weekly binned light curve, we obtained a 
light curve with 30 d integrations using the unbinned likelihood 
analysis to account for the effect of sparse sampling of certain 
sources in our radio sample during the correlation. The rest of the 
procedure to obtain the monthly binned light curve is similar to 
those of weekly binned described above. 

3.2 Metsahovi 

The 37 GHz observations were obtained with the 13.7 m diameter 
Metsahovi radio telescope, which is a radome enclosed paraboloid 
antenna situated in Finland. The measurements were made with a 
1 GHz-band dual beam receiver centered at 36.8 GHz. The obser¬ 
vations are ON-ON observations, alternating the source and the 
sky in each feed horn. A typical integration time to obtain one flux 
density data point is between 1200 and 1400 s. The detection limit 
of our telescope at 37 GHz is on the order of 0.2 Jy under optimal 
conditions. Data points with a signal-to-noise ratio < 4 are handled 
as non-detections and discarded from the analysis. 

The flux density scale is set by observations of the HII re¬ 
gion DR2L with a known flux density of 17.9 Jy at 37 GHz 
(Teriisranta et al. 1998). Sources NGC 7027, 3C 274 and 3C 84 are 
used as secondary calibrators. A detailed description of the data 
reduction and analysis is given in Terasranta et al. (1998). The er¬ 
ror estimate in the flux density includes the contribution from the 
measurement rms and the uncertainty of the absolute calibration. 


4 THE DISCRETE CORRELATION FUNCTION 

The cross-correlation function is commonly employed in the study 
of AGN to probe the continuum emission mechanism by corre¬ 
lating multiwavelength light curves and to seek correlations be¬ 
tween the variability and other AGN properties. The classical cor¬ 
relation function by Oppenheim & Schafer (1975) can be used for 
such purposes, if the data are evenly sampled. To deal with the un¬ 
even sampling, there are three variants of the classical function - 
Interpolated cross-correlation function (Gaskell & Peterson 1987), 
discrete correlation function (DCF; Edelson & Krolik 1988) and 
^-transformed DCF (Alexander 2013). Due to unevenly sampled 
7 -ray and radio light curves of sources in this work we resorted to 
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discrete correlation function. The DCF is estimated from the rela¬ 
tion, 

DCF., = (j) 

where ai, bj are the observed fluxes at times ti and tj and a, b, a a 
and CTf, are the means and standard deviation of the entire light 
curves. At this step, DCF;, are binned by their associated time lag, 
Tij = ti — tj into equal width bins. The average of the bins yields 
the DCF(t). 

The normalization (mean and standard deviation of respec¬ 
tive light curves) used in the estimation of the DCF(r) assumes 
that the light curves are statistically stationary. Under this assump¬ 
tion, the DCF sometimes can exceed unity (DCF(r) > 1) mak¬ 
ing the interpretation difficult. A workaround for this problem is 
to estimate the mean and standard deviations only from the points 
that overlap at a given time lag bin (c.f.. White & Peterson 1994; 
Welsh 1999). Thus, with the local normalization the resultant DCF 
is bound to the [—1, -fl] interval. A positive DCF(r) implies a cor¬ 
related variability and an anti-correlation when its negative. The 
uncertainties of the DCF were estimated by a model-independent 
Monte Carlo method (Peterson et al. 1998) that accounts for the 
effects of measurement noise and data sampling. The simulation 
consists of a bootstrap selection of a subsample of data points from 
each light curve to which Gaussian noise with a standard deviation 
matching the observational error bars are added. The time lags and 
68 % fiducial interval are estimated following the maximum likeli¬ 
hood approach by Alexander (2013). According to this method, the 
maximum likelihood estimate coincides with the peak of the cross¬ 
correlation. The fiducial interval is then estimated by interpolating 
between the points of the likelihood function (see Alexander 2013). 

Works such as, Welsh (1999) have shown that removing a lin¬ 
ear trend from the light curve prior to estimating the DCF to im¬ 
prove the estimation of time lag. However, this was found not to 
be true in the case of unevenly sampled light curves (Peterson et al. 
2004). 


5 SIGNIFICANCE OF THE CORRELATIONS 

The stochastic nature of the variability, data sampling and measure¬ 
ment noise are some of the properties that affect the result of the 
correlation. Also, the frequent appearance of flares means that high 
correlation coefficients between any two energy bands are to be 
expected even in the absence of any physical relation between the 
processes responsible for their production. Hence, we investigated 
the statistical significance of the DCF by performing Monte Carlo 
simulations, in turn, estimating the probability that the observed 
correlation is primarily limited by chance correlations. The widely 
employed procedure to estimate the significance is from the cross¬ 
correlation of simulated light curves with power law power spectral 
densities. Recently, two additional methods have been proposed 
(Fuhrmann et al. 2014) to quantify the significance of the cross¬ 
correlation using mixed source correlations and stacking analysis. 

We adopted all the three methods in estimating the signifi¬ 
cance and have discussed the methods in detail below. Every corre¬ 
lation is locally normalized as some works have shown the method 
to be efficient also with simulated light curves (e.g., Welsh 1999; 
Max-Moerbeck et al. 2014b). 


5.1 Correlation of simulated light curves 

The AGN light curves, in general, can be modelled by red- 
noise power spectra showing variability at all time scales, 
e.g., Hufnagel & Bregman (1992) in the radio and optical, 
Lawrence & Papadakis (1993) in the X-rays, and Abdo et al. 
(2010) in 7 rays. Hence to construct the significance level of the 
observed correlations, we correlate simulated light curves for both 
radio and 7 rays assuming a simple power law model of the form: 
PSDoc /-“. 

5. 1.1 Estimation of the Power Spectral Density 

The underlying PSD of the light curve is estimated by fitting a sim¬ 
ple power law to the periodogram of the observed light curve or by 
fitting a straight-line to the periodogram in logarithmic scale. Due 
to the complications arising from red-noise leak and aliasing, we 
estimated the PSD using a variant of PSRESP method proposed by 
Uttley et al. (2002). This method involves the simulation of light 
curves for a set of model parameters, estimate the periodogram for 
each simulated light curve and average them, and determine the 
best-fitting PSD from a goodness of fit. We discuss the method in 
detail in the following sub-sections. 

Estimation of the Periodogram. To estimate the periodogram, every 
light curve for k = 1,2,... ,N) is initially binned in time 

intervals ATbin, by taking the weighted mean of all points within 
each bin. Missing data points in the binned light curve were linearly 
interpolated before convolving with an Hanning window function 
to reduce red-noise leakage (Max-Moerbeck et al. 2014b). 

For an evenly sampled light curve of length N, the peri¬ 
odogram is given by the modulus-squared of the discrete Fourier 
transform (DFT; Press et al. 1992) of the data. The periodogram, 
P{fj), at a given Fourier frequency fj is then given by, 

Pifi) = {Re[DFT{j)f + Im[DFTij)f} , (2) 

where NAThin/ is the fractional rms normalization. With this 
normalization, the square root of the integral of the underlying PSD 
between two frequencies fi and /2 yields the contribution to the 
fractional rms squared variance (i.e., ffd) due to variations from 
the corresponding time-scales (Miyamoto et al. 1991; van der Klis 
1997). Thus, integration between /i and /Nyq (even) or /(„_i )/2 
(odd) yields the total rms squared variability. 

The resulting periodogram is then logarithmically transformed 
and binned with 10 points as discussed by Papadakis & Lawrence 
(1993). This approach reduces the scatter in the periodogram. 
But taking logarithm also introduces a bias to the power of 
the periodogram which is a constant and can be removed by 
adding 0.25068 to obtain the final binned logarithmic periodogram 
(Vaughan 2005). The binned logarithmic periodogram are normally 
distributed within each geometric mean frequency bin. 

Simulating Light Curves. To simulate a light curve, we specify a 
power spectral model which we wish to test against the data (power 
law in our case). The normalization of the model power spectrum 
is a multiplicative factor which is carried through any convolution 
with the window function (i.e., only the power spectral shape is 
distorted by sampling). 

Due to the finite length of the observed light curve, power 

® fj = j! (AATbin). The zero Fourier frequency component, /q = 0, cor¬ 
responds to the sum of the light curve. For even N, f 1 , 1/2 = U(2 ATbin) is 
the Nyquist frequency, /Nyq- 
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from longer than observed time-scales leaks into the shorter time- 
scales and distorts the observed PSD. This effect, called red-noise 
leak, was accounted for by simulating light curves 100 times longer 
than the observed light curve. 

On the other hand, uneven sampling of the light curve causes 
variations on time-scales down to the resolution of the observed 
light curve by aliasing. The power above /Nyq which is aliased 
to frequencies below /Nyq makes the observed periodogram dis¬ 
torted. This is taken into account by limiting the resolution of the 
simulated light curves to be O.lATsamp. 

By considering the above complications, we simulate N light 
curves assuming a power law using the method proposed by 
Emmanoulopoulos et al. (2013) (hereafter EMP13). The widely 
used algorithm of Timmer & Koenig (1995) (hereafter TK95) for 
simulating light curves is appropriate for the production of Gaus¬ 
sian light curves. Since the distribution of most of the light curves 
in our sample are far from being a Gaussian, light curves simulated 
using the method by TK95 may not be appropriate for the establish¬ 
ment of confidence intervals for PSD and cross-correlation studies. 
The method by EMP13 involves the combination of the routine by 
TK95 and the iterative amplitude adjusted Eourier transform algo¬ 
rithm by Schreiber & Schmitz (1996), producing light curves pos¬ 
sessing exactly the PSD and the probability density function (PDE) 
as the observed light curve. Eor more information on the method 
the reader is referred to Emmanoulopoulos et al. (2013). 

We simulated 1000 continuous light curves using the method 
by EMP13. Instead of generating random numbers from the best¬ 
fitting PDE of the data, which is used during the amplitude adjust¬ 
ment stage in EMP13, we considered the use of cumulative dis¬ 
tribution function (CDF). This was due to the multi-modal nature 
of the observed data making it very difficult to fit a function to 
the PDF. It is not quite so simple to estimate a PDF. If one uses 
a histogram one needs to choose the bin width and the starting 
point for the first bin. If one use the kernel density estimation one 
needs to choose the kernel shape and bandwidth. Hence, to circum¬ 
vent all the problems of obtaining a PDF, we generated the random 
numbers from the CDF of the data. The CDF has a simple non- 
parametric estimator that needs no choices to be made - the empir¬ 
ical distribution function. We also confirmed that the distribution 
of the generated random numbers agree with the observed data by 
using a two-sample Kolmogorov-Smimov (KS) test. 

The simulated light curve is finally normalized to the mean 
and variance of the observed light curve. 

Goodness of fit. Once a continuous light curve is simulated, it is 
resampled to the sampling pattern as the observed light curve. At 
this stage, Gaussian noise with zero mean and variance matching 
those of the observations are added to the resampled light curve in 
the case of radio light curves, while for the 7 rays Poisson noise is 
added according to the equation. 


LCsim (ft 


Pois\pi = 

At 


for i = 1,..., A, (3) 


where Pois[tL\ corresponds to the Poisson random number with 
a mean value of LCsiin{ti)At. The resampled light curve is then 
rebinned, and empty bins interpolated in the same manner as for 
the observed light curve. The periodogram of the resulting light 
curve is then obtained using equation 2. We thus obtained 1000 
periodograms for every simulated light curve which was binned 
logarithmically as was implemented for the observed periodogram. 
The binned logarithm periodograms are then averaged to obtain the 
model average periodogram, Psim(/), and the rms spread about the 


mean is also calculated and taken as the error in the power at each 
frequency, APaimif). 

Having obtained the model periodogram, we now estimate a 
statistic defined as Xdist (Uttley et al. 2002) calculated from the 
model and observed periodogram, Poba(/) as follows. 


2 ^ [Psimjf) - Pobs(/)]^ 

APaim(/)2 


(4) 


Next, we determine the Xdist between the model and every simu¬ 
lated periodogram by replacing the Pobs(/) with Psim,i(/) (where 
i = 1,2,... ,1000) in equation 4. The goodness of fit is then given by 
the percentile of the simulated Xdist distribution that exceeds the 
observed Xdist estimate. 

This whole approach from simulating light curves to estimat¬ 
ing the goodness of fit was tested for a range of PSD values. The 
PSD value with the highest goodness of fit corresponds to the best¬ 
fitting PSD. Using this PSD value we simulate 1000 light curves 
using the same approach as discussed above, obtain their binned 
logarithm periodogram and estimate the PSD for each periodogram 
by linear least squares method. The distribution of PSD obtained 
was used to get the la confidence interval for the best-fitting PSD. 

Thus having obtained the best-fitting PSD of all the light 
curves in our sample, we simulate 1000 light curves using the al¬ 
gorithm proposed by EMP13 for both radio and 7 rays. To every 
simulated light curve assuming the power law properties, we added 
Gaussian noise matching those of the observations for the case of 
radio, while for the 7 rays Poisson noise was added as shown in 
equation 3. The cross-correlation function of the simulated light 
curves were then estimated using the same approach as for the ob¬ 
served light curve. 

This simulation is followed by estimating the 68.27% (Icr), 
95.45% {2a) and 99.73% (3 ct) significance levels from the distri¬ 
bution of the cross-correlation coefficients for each time lag bin. 


5.2 Mixed source correlations 

The significance levels of the cross-correlation under this context 
is estimated by correlating all sources with possible combinations 
from our sample excluding the source under study. That is, having 
55 sources in the source list, we correlate every source in the radio 
with all the 54 7 -ray light curves, in turn, yielding 2916 correla¬ 
tions. The correlations are done just as for the real data. The 1, 2 
and 3a significance levels are estimated similar to that discussed 
in the above section. This approach of using every source in the 
sample instead of simulating light curves is under the assumption 
that the flares exhibited by the source at different wavelengths are 
physically unrelated. It is similar to the method where light curves 
are simulated over a range of PSD under the assumption that all the 
sources exhibit similar variability properties, i.e., characterized by 
red-noise PSDs (e.g., Agudo et al. 2011a,b; Schinzel et al. 2011). 
This also assumes that the light curves are sampled in the same 
way, which is strictly not true in the case of our radio observations. 


5.3 Stacking the correlations 

Following Fuhrmann et al. (2014), to improve the sensitivity for 
the detection of correlations we consider stacking or averaging the 
correlations obtained from the whole source sample. We also at¬ 
tempted to stack the correlations obtained from the light curves 
which were initially normalized by dividing with the mean flux 
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density. For further details on this method, see Fuhrmann et al. 
(2014). 

6 RESULTS 

6.1 Power Spectral Density 

The PSDs characterize the variability time-scale of a source. For 
estimating the PSDs of the radio light curves, we considered all 
the data available until 2013 August 3 to account for the long¬ 
term variability exhibited by the source. In most of the cases, we 
used over 20 yr of data from the Metsahovi monitoring programme 
(Terasranta et al. 2005). We estimated the PSD slopes of every 
source at both frequencies following the Monte Carlo simulation 
as discussed in Section 5.1.1. From the simulation we obtained the 
acceptance probability (p) which allowed us to set a constraint such 
that only the PSD slopes of sources with p > 0.05 are considered 
to be robust estimates. The PSD and its acceptance probability at 

7 rays and radio for 3C 454.3 are shown in Fig. 1. In radio we 
obtained the PSD slopes for 51 sources, while in 7 rays only 48 
sources had acceptable estimates. The PSD slopes of other sources 
are taken from the average of the PSD slopes, which in radio for 
FSRQs and BL Lacs are 2.2 and 2.0, respectively. In 7 rays the 
mean PSD slope for FSRQs and BL Lacs are 1.3 and 1.1, respec¬ 
tively. This indicates that the long-term variability dominates the 
radio light curves, while in the 7 rays the short-term variations are 
more pronounced. The PSD slopes for individual sources are shown 
in Table 1 . The distribution of the PSD slopes in radio and in 7 rays 
are shown in Fig. 2. 

Based on a comparison of our PSD estimates with earlier find¬ 
ings from the literature, we find: 

(i) 3C 273: In the 7 rays, Nakagawa & Mori (2013) and 
Sobolewska et al. (2014) reported 1.30 [1.04, 1.56] and 0.84 [0.75, 
0.95], respectively, which are consistent given the errors in this 
work. 

(ii) 3C 279: In the radio, Chatterjee et al. (2008) reported a PSD 
slope of 2.3. Considering our higher frequency radio data and the 
fact that the short time-scale variability increases with frequency, 
our estimate of 2.00 [1.93, 2.11] is consistent with theirs. 

(iii) BL Lac: Sobolewska et al. (2014) in the 7 rays obtained 
0.93 [0.79, 1.11], consistent with 1.14 [1.10, 1.22] reported in this 
work. 

(iv) 3C 454.3: In the 7 rays, Nakagawa & Mori (2013) obtained 
1.49 [1.33, 1.65] in good agreement within the Icr confidence in¬ 
terval reported in this work. 

(v) Abdo et al. (2010) using first 11 months of FermifLfd data 
reported an averaged PSD slope for 22 FSRQs and 6 BL Lacs to 
be 1.5 and 1.7, respectively. These are higher than our estimates 
obtained from 35 FSRQs and 12 BL Lacs. This difference could 
possibly be due to the longer time range of FermifLPtT data con¬ 
sidered in this work. However, the PSD slope of 1.6 [1.4, 1.8] for 
3C 279 reported by the authors is consistent with our result within 
the quoted Icr confidence interval. 


6.2 Cross-Correlation significance in individual sources 

The cross-correlation of weekly/monthly binned 7 -ray and radio 
light curves were performed as discussed in Section 4. Unlike 
the 7 -ray light curves, the radio light curves of some sources 
show an increasing or decreasing long-term trend which results in 


larger cross-correlation coefficients making the interpretation diffi¬ 
cult. Six sources in our sample - 3C 84, 0458—020, 0605—085, 
0917-1-449, 3C 446 and BL Lac - were found to exhibit such a 
trend. These sources were linearly detrended prior to the correla¬ 
tion analysis. When investigating the time lags, we only considered 
lags up to half of the duration of the shortest light curve, in order 
to avoid spurious correlations. The significance of the DCF peak 
is estimated by simulating light curves of known PSDs (hereafter 
method 1; Section 5.1) and also from mixed source correlations 
(hereafter method 2; Section 5.2). 

Cross-correlating weekly binned light curves and using 
method 1 to estimate the significance, we find 23 and 10 significant 
correlations at the 2a and 3a levels, respectively. With method 2 we 
find 38 and 6 sources significant above the 2a and 3cr levels, respec¬ 
tively. From the cross-correlation of monthly binned light curves, 
we obtained 23 and 13 sources significant at 2a and 3cr levels, re¬ 
spectively, using method 1. The number of significant sources at 2a 
and 3cr levels were found to be 30 and 6 , respectively, using method 
2 . 

Following Max-Moerbeck et al. (2014a), we wish to set a 
threshold, which would produce only one spurious correlation in 
our sample based on chance probability. For our sample (55 sources 
in total), this is 98.18% (2.36cr). We found 20 sources from the 
correlation of weekly binned and 23 sources from monthly binned 
7 -ray with radio light curves significant at 2.36cr level. The 
time lag and the DCF peak for sources with correlation significance 
>2.36(7 level are given in Table 2. A positive time lag denotes that 
the radio lags the 7 rays. The distribution of the lags of sources re¬ 
ported in Table 2, using both weekly and monthly light curves, are 
shown in Fig. 3. 

We find that sources with higher significance level have the 
sampling of the 7 -ray light curve closer to the sampling of the ra¬ 
dio light curve, justifying the use of two binning schemes (weekly 
and monthly binned). However, when the significance level was the 
same we found no major difference in the time lags from the two 
binned cases. In few cases, were such a difference exists we sug¬ 
gest to consider the result obtained from the weekly binned 7 -ray 
light curves. 

6.3 Stacking the correlations 

The stacked DCFs were obtained following the method described 
in Section 5.3. We got identical results from both the methods dis¬ 
cussed. Hence, we adhered to the first method of simply averaging 
the correlations. We stacked the DCFs for the whole sample and 
also for FSRQs and BL Lacs, respectively. The significance level 
for the stacked DCFs are estimated using the mixed source corre¬ 
lations. The stacked DCFs for the whole sample in observer’s and 
source’s frame, by scaling the time lags with a factor of l/(l-|- 2 :), 
are shown in Fig. 4, while the stacked DCFs for the sub-samples are 
shown in Fig. 5. The results obtained for the whole sample in the 
observer’s frame are 80 and 120 d for weekly and monthly binned 
light curves while for those in source frame are 47 and 70 d, re¬ 
spectively (see Table 3). 

The DCFs obtained for the whole sample and for the sub¬ 
samples are significant at > 3a level. However, due to the fewer 
BL Lacs in our sample, the errors of their stacked DCFs are higher. 
Owing to the broad DCF peak, we cannot distinguish between the 
time lags obtained between the weekly and monthly-binning. 

To test for the possible bias on the significance of the stacked 
DCFs that might have been introduced by the sample selection, we 
removed all the sources with significance ^ 2.36a. The peak of 
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3C 454.3: PSD in radio 


3C 454.3: PSD in 7 -ray 





Figure 1. Power Spectral Density of the blazar 3C454.3 in 7 rays (Top left) and in radio (Top right) and their corresponding acceptance probabilities in the 
bottom panel. The PSD slope is shown in the figure along with its Icr confidence interval in square brackets. 



Figure 2. Distribution of PSD slopes in radio and in 7 rays. 
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Table 1. Source sample and PSD results. 


Source 

(1) 

2FGL Name 
(2) 

Optical class 
(3) 

z 

(4) 

PSDradio 

(5) 


PSD.^ 

(6) 

0059+581 

J0102.7+5827 

FSRQ 

0.644 

1.53 

[1.51, 1.61] 

1.75 

[1.70, 1.85] 

0106+013 

J0108.6+0135 

FSRQ 

2.099 

2.08 

[2.05, 2.13] 

1.11 

[1.03, 1.17] 

0133+476 

J0136.9+4751 

FSRQ 

0.859 

2.31 

[2.29, 2.39] 

1.05 

[1.00, 1.12] 

0212+735 

J0217.7+7353 

FSRQ 

2.367 

1.97 

[1.84,2.15] 

0.48 

[0.43, 0.55] 

0215+015 

J0217.9+0143 

FSRQ 

1.721 


- 

0.57 

[0.51,0.63] 

0234+285 

J0237.8+2846 

FSRQ 

1.206 

2.52 

[2.45, 2.70] 

1.65 

[1.60, 1.76] 

0235+164 

J0238.7+1637 

BL Lac 

0.940 

2.23 

[2.17, 2.43] 

1.15 

[1.06, 1.31] 

0300+470 

J0303.5+4713 

BL Lac 

- 


- 

0.63 

[0.55, 0.68] 

3C 84 

J0319.8+4130 

GAL 

0.017 

1.54 

[1.48, 1.58] 

1.13 

[1.09, 1.21] 

CTA 026 

J0339.4-0144 

FSRQ 

0.852 

2.28 

[2.22, 2.44] 

1.42 

[1.39, 1.51] 

0420-014 

J0423.2-0120 

FSRQ 

0.916 

1.95 

[1.79, 2.00] 

1.81 

[1.74, 1.86] 

0440-003 

J0442.7-0017 

FSRQ 

0.844 

2.10 

[2.01,2.34] 

1.67 

[1.63, 1.75] 

0458-020 

J0501.2-0155 

FSRQ 

2.291 

1.72 

[1.62, 1.93] 

0.84 

[0.79, 0.90] 

0528+134 

J0530.8+1333 

FSRQ 

2.070 

2.21 

[2.13, 2.33] 

0.82 

[0.78, 0.89] 

0605-085 

J0608.0-0836 

FSRQ 

0.870 

2.23 

[1.99, 2.31] 

1.07 

[1.02, 1.14] 

0716+714 

J0721.9+7120 

BL Lac 

0.310 

2.00 

[1.92, 2.15] 

0.99 

[0.90, 1.06] 

0736+017 

J0739.2+0138 

FSRQ 

0.189 

1.79 

[1.75, 1.89] 

1.12 

[1.07, 1.20] 

0754+100 

J0757.1+0957 

BL Lac 

0.266 

1.90 

[1.84, 2.06] 


- 

0805-077 

J0808.2-0750 

FSRQ 

1.837 


- 

1.00 

[0.96, 1.11] 

0814+425 

J0818.2+4223 

BL Lac 

- 

2.54 

[2.43, 2.68] 

0.84 

[0.78, 0.90] 

OJ248 

J0830.5+2407 

FSRQ 

0.942 

2.50 

[2.35, 2.65] 

1.52 

[1.49, 1.64] 

0836+710 

J0841.6+7052 

FSRQ 

2.218 

1.76 

[1.71, 1.87] 

1.82 

[1.77, 1.91] 

OJ 287 

J0854.8+2005 

BL Lac 

0.306 

2.14 

[2.04, 2.29] 

1.12 

[1.06, 1.23] 

0917+449 

J0920.9+4441 

FSRQ 

2.188 

2.69 

[2.50, 2.97] 

1.66 

[1.60, 1.77] 

0954+658 

J0958.6+6533 

BL Lac 

0.367 

1.73 

[1.64, 1.87] 


- 

1055+018 

J1058.4+0133 

FSRQ 

0.888 

2.19 

[2.16, 2.28] 

1.59 

[1.56, 1.71] 

1150+497 

J1153.2+4935 

FSRQ 

0.333 

2.56 

[2.42, 2.87] 

1.84 

[1.79, 1.95] 

1156+295 

J1159.5+2914 

FSRQ 

0.725 

2.34 

[2.31, 2.42] 

1.37 

[1.33, 1.47] 

1222+216 

J 1224.9+2122 

FSRQ 

0.434 

1.60 

[1.49, 1.68] 

1.49 

[1.43, 1.58] 

3C 273 

J1229.1+0202 

FSRQ 

0.158 

2.27 

[2.18, 2.37] 

1.18 

[1.08, 1.22] 

3C 279 

J1256.1-0547 

FSRQ 

0.536 

2.00 

[1.93, 2.11] 

1.42 

[1.36, 1.53] 

1308+326 

J1310.6+3222 

FSRQ 

0.997 

2.08 

[2.03, 2.20] 

1.27 

[1.23, 1.36] 

1334-127 

J1337.7-1257 

FSRQ 

0.539 

2.33 

[2.24, 2.50] 


- 

1502+106 

J1504.3+1029 

FSRQ 

1.839 

2.28 

[2.11, 2.41] 

1.15 

[1.11, 1.28] 

1510-089 

J1512.8-0906 

FSRQ 

0.360 

2.39 

[2.33, 2.54] 

1.38 

[1.30, 1.44] 

1546+027 

J 1549.5+0237 

FSRQ 

0.414 

2.84 

[2.76, 3.24] 


- 

1633+382 

J1635.2+3810 

FSRQ 

1.813 

2.25 

[2.16, 2.39] 

1.42 

[1.37, 1.53] 

1638+398 

J 1640.7+3945 

FSRQ 

1.666 

2.24 

[1.97, 2.33] 

2.14 

[2.11,2.26] 

3C 345 

J 1642.9+3949 

FSRQ 

0.593 

1.81 

[1.78, 1.84] 


- 

Mai-k 501 

J1653.9+3945 

BL Lac 

0.033 

1.44 

[1.41, 1.49] 

1.76 

[1.71, 1.82] 

1730-130 

J1733.1-1307 

FSRQ 

0.902 

1.77 

[1.75, 1.84] 

1.12 

[1.07, 1.23] 

1739+522 

J 1740.2+5212 

FSRQ 

1.379 


- 

0.83 

[0.79, 0.90] 

1749+096 

J1751.5+0938 

BL Lac 

0.322 

2.27 

[2.25, 2.33] 

0.60 

[0.54, 0.66] 

1803+784 

J1800.5+7829 

BL Lac 

0.680 

1.89 

[1.84, 1.97] 

1.01 

[0.97, 1.08] 

3C 371 

J 1806.7+6948 

BL Lac 

0.051 

1.88 

[1.83, 1.99] 

1.20 

[1.16, 1.27] 

4C 56.27 

J 1824.0+5650 

BL Lac 

0.664 

2.13 

[2.07, 2.29] 

1.34 

[1.31, 1.43] 

1828+487 

J 1829.7+4846 

FSRQ 

0.692 

2.40 

[2.31, 2.50] 


- 

2022-077 

J2025.6-0736 

FSRQ 

1.388 

2.35 

[2.23, 2.48] 

1.36 

[1.32, 1.45] 

BL Lac 

J2202.8+4216 

BL Lac 

0.068 

2.02 

[1.96, 2.04] 

1.14 

[1.10, 1.22] 

2201+171 

J2203.4+1726 

FSRQ 

1.076 

2.50 

[2.24, 2.56] 

0.90 

[0.85, 0.96] 

3C446 

J2225.6-0454 

FSRQ 

1.404 

2.17 

[2.14, 2.23] 


- 

2227-088 

J2229.7-0832 

FSRQ 

1.559 

2.16 

[2.10, 2.39] 

1.24 

[1.19, 1.33] 

2230+114 

J2232.4+1143 

FSRQ 

1.037 

2.31 

[2.22, 2.51] 

0.97 

[0.88, 1.03] 

2234+282 

J2236.4+2828 

BL Lac 

0.795 

1.88 

[1.74, 1.95] 

1.69 

[1.64, 1.76] 

3C 454.3 

J2253.9+1609 

FSRQ 

0.859 

2.31 

[2.12, 2.44] 

1.65 

[1.60, 1.72] 


Columns are as follows: (1) Source name; (2) 2FGL name; (3) optical classification; (4) redshift 
(Nolan et al. 2012); (5){6) the first value is the best-fitting PSD obtained from Monte Carlo simulations 
at radio and 7 -ray frequencies along with the 68.27% confidence intervals in square brackets. 
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Figure 3. Distribution of Lags in weekly (Left) and monthly binned (Right). 


Table 2. Time lag and correlation of individual sources with significance >2a. The significance level (tr) of the DCFpeak 2 nd the distance travelled by the 
emission region (d-^.radio) in parsecs are also shown. Time lags with negative sign denotes that the radio leads the 7 -ray, and vice versa for time lags with 
positive sign. 


Weekly Monthly 


Source 

Lag 

(d) 

DCF 

a 

^7,radio 

(pc) 

Lag 

(d) 

DCF 

a 

^7,radio 

(pc) 

0059-H581 

160 [152, 172] 

0.766 

>3 

11.3 [10.7, 12.2] 

140 [130, 176] 

0.842 

>3 

9.9 [9.2, 12.5] 

0106-H013 

- 

- 

- 

- 

350 [205, 446] 

0.640 

>2.36 

44.3 [26.0, 56.6] 

0133-H476 

150 [142, 163] 

0.568 

>2.36 

22.6 [21.4, 24.5] 

150 [120, 163] 

0.662 

>2.36 

22.6 [18.1, 24.5] 

0215-H015 

- 

- 

- 

- 

600 [573, 720] 

0.652 

>2.36 

- 

0234-H285 

60 [50, 96] 

0.614 

>2.36 

8.5 [7, 13.6] 

90 [60, 112] 

0.651 

>2.36 

12.7 [8.5, 15.8] 

0235-H164 

30 [16, 36] 

0.928 

>3 

14.3 [7.6, 17.2] 

20 [11, 29] 

0.959 

>3 

9.5 [5.2, 13.8] 

3C 84 

500 [493, 508] 

0.491 

>2.36 

- 

- 

- 

- 

- 

0420-014 

- 

- 

- 

- 

20 [8, 30] 

0.615 

>2.36 

1.9 [0.7, 2.8] 

0440-003 

-540 [-591, -505] 

0.618 

>2.36 

-20.4 [-22.3, -19.1] 

- 

- 

- 

- 

0805-077 

120 [97, 156] 

0.590 

>2.36 

- 

100 [78, 133] 

0.753 

>3 

- 

0814-H425 

210 [128, 230] 

0.524 

>2.36 

- 

200 [151,234] 

0.703 

>2.36 

- 

OJ248 

-720 [-738, -702] 

0.829 

>3 

-88.0 [-90.2, -85.8] 

-720 [-736, -707] 

0.916 

>3 

-88.0 [-89.9, -86.4] 

1156-H295 

- 

- 

- 

- 

25 [11, 92] 

0.617 

>2.36 

8.6 [3.8, 31.8] 

1222-H216 

280 [261,287] 

0.601 

>2.36 

60.4 [56.3,61.9] 

260 [247, 290] 

0.703 

>2.36 

56.1 [53.3,62.6] 

3C 273 

160 [151, 165] 

0.625 

>2.36 

29.6 [27.9, 30.6] 

160 [149, 169] 

0.770 

>2.36 

29.6 [27.6, 31.3] 

1308-H326 

- 

- 

- 

- 

425 [412, 498] 

0.694 

>2.36 

89.0 [86.2, 104.3] 

1502-H106 

30 [23, 88] 

0.872 

>3 

2 [1.5, 5.9] 

50 [38, 63] 

0.940 

>3 

3.3 [2.5, 4.2] 

1633-H382 

-15 [-21, -4] 

0.686 

>3 

-2.9 [-4.1, -0.8] 

0 [-27, 11] 

0.701 

>2.36 

0 [-5.3, 2.2] 

3C 345 

-40 [-61, -11] 

0.561 

>2.36 

-4.6 [-7.0, -1.3] 

30 [-3, 68] 

0.648 

>2.36 

3.4 [0.3, 7.8] 

3C 345“ 

80 [47, 101] 

0.551 

>2.36 

9.2 [5.4, 11.6] 

- 

- 

- 

- 

1730-130 

100 [75, 110] 

0.666 

>2.36 

19.1 [14.3,21.0] 

40 [31, 52] 

0.798 

>3 

7.6 [5.9, 9.9] 

1749-H096 

135 [120, 144] 

0.513 

>2.36 

8.8 [7.8, 9.4] 

25 [14, 42] 

0.568 

>2.36 

1.6 [0.9, 2.7] 

2022-077 

-850 [-861, -838] 

0.678 

>2.36 

- 

- 

- 

- 

- 

BL Lac 

- 

- 

- 

- 

-620 [-641, -600] 

0.638 

>2.36 

-35.3 [-36.6, -34.2] 

2201-H171 

560 [522, 624] 

0.609 

>2.36 

- 

690 [541,703] 

0.744 

>2.36 

- 

2230-H114 

- 

- 

- 

- 

550 [529, 576] 

0.662 

>2.36 

36.1 [34.7, 37.7] 

3C 454.3 

40 [32, 61] 

0.621 

>2.36 

11.7 [9.3, 17.8] 

40 [30, 59] 

0.732 

>2.36 

11.7 [8.7, 17.2] 


“ second DCF peak for 3C 345 with less significance than the first although being above the 2.36 (t level 


the resultant stacked DCFs was lower by a factor of ~ 1.2 when 
compared to the overall stacked DCFs but was still significant at 
> 3a level. This remains the same in the case of FSRQs while for 
the BL Lacs the stacked DCFs were significant only at 90%. Thus, 
the significance of the stacked DCFs is not affected by sources with 
significant correlation. 


7 DISCUSSION 

7.1 Correlations and comparison with earlier results 

Observations have shown that 7 -ray loud AGNs are clearly asso¬ 
ciated with compact, flat radio spectrum sources (Ackermann et al. 
2011). The emission and variability from the relativistic jet in the 
radio and the 7 rays may both be coupled to the disturbances in 
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stacked DCF: weekly-binned (observer’s frame) Stacked DCF: monthly-binned (observer’s frame) 




Stacked DCF: weekly-binned (source’s frame) Stacked DCF: monthly-binned (sonrce’s frame) 




Figure 4. Stacked DCFs obtained from correlation of radio and weekly binned (Left column) and likewise from radio and monthly binned (Right column) 
7 -ray light curves. Stacked DCFs for the whole sample shown in observer’s frame (Top) and source frame (Bottom). The red, green and blue dotted lines 
correspond to Icr, 2a and 3cr significance levels, respectively. 


Table 3. Results from stacking the correlation both in observers and source frame. Columns are similar to Table 2. 





Weekly 



Monthly 


Sample Type 

Frame 

Lag 

DCF 

^" 7 , radio 

Lag 

DCF 

^ 7 ,radio 



(d) 


(pc) 

(d) 


(pc) 

Whole sample 

Observer 

80 [7, 206] 

0.200 

7.0 [0.6, 18.1] 

120 [41,216] 

0.250 

10.6 [3.6, 19.0] 

Source 

47 [4, 125] 

0.219 

4 [0.3, 11] 

70 [23, 134] 

0.280 

6 [2, 12] 

FSRQs 

Observer 

80 [-1,222] 

0.230 

7.0 [-0.1, 19.6] 

90 [13, 212] 

0.306 

7.9 [1.1, 18.6] 

Source 

47 [2, 141] 

0.229 

4 [0.1, 12] 

53 [7, 128] 

0.306 

4.7 [0.6, 11.3] 

BL Lacs 

Observer 

120 [93, 274] 

0.200 

10.6 [8.2,24.1] 

120 [29, 157] 

0.220 

10.5 [2.5, 13.8] 

Source 

62 [-22, 83] 

0.199 

5.5 [-2, 7.3] 

62 [-22, 106] 

0.223 

5.5 [-2, 9.3] 


the central engine. The flares seen in the radio light curves are 
physically linked to the ejections of superluminal radio components 
(e.g., Savolainen et al. 2002; Fromm et al. 2013). Therefore, cross¬ 
correlation and time lags between radio and 7 -ray light curves can 


be used to place constrains on the location of the 7 -ray flares as 
they cannot be spatially resolved by the existing instruments. 

For many sources in this work we have obtained significant 
correlations as shown in Section 6 . From the correlation results ob¬ 
tained using weekly binned 7 -ray light curves, we found 3 sources 


MNRAS 000, 1-15 (2015) 












































































10 


Ramakrishnan et al. 


stacked DCF: FSRQ (weekly-binned) Stacked DCF: FSRQ (monthly-binned) 




Stacked DCF: BL Lac (weekly-binned) Stacked DCF: BL Lac (monthly-binned) 




Figure 5. Stacked DCFs for the sub-sample comprising of FSRQs (Top) and for sub-sample comprising of BL Lacs (Bottom) are shown. In left column shown 
are the correlation of radio and weekly binned and those from radio and monthly binned 7 -ray are shown in right column. 


(0440—003, OJ 248 and 2022—077) displaying a negative time lag 
(radio leading) with At >1.5 yr. Likewise, using monthly binned 
7 -ray light curves, two sources (OJ 248 and BL Lac) showed a 
negative time lag with At > 1.5 yr. In all of these cases, except 
for BL Lac, there is a 7 -ray flare at the very beginning or end of 
the 5 yr period and it is difficult to judge the reliability of the cor¬ 
relation without complete sampling of the flares. In BL Lac the 
monthly binned 7 -ray light curve averages the variations too much 
and the weekly binning (where the correlation is not significant) is 
closer to the sampling of the radio light curve. Therefore, these re¬ 
sults are not considered for further interpretation, in turn reducing 
the number of significant cases to 16 and 21 sources for weekly and 
monthly binned, respectively. For 15 sources in our sample, this is 
the first time a significant correlation is reported in the literature. 
Our results agree very well with those quoted in Table 4 from ear¬ 
lier works. 


Table 4. Time lags obtained from radio/ 7 -ray correlation in earlier works. 


Source 

Lag 

(d) 

Frequency 

(GHz) 

Reference 

0234-H285 

40 [30, 50] 

86 

Fuhrmann et al. (2014) 

0235-H164 

0-50 

15,43 

Agudo et al. (201 lb) 


-4 [-14, 6 ] 

86 

Fuhrmann et al. (2014) 

3C 273 

120-170 

43 

Beaklini & Abraham (2014) 

1502-H106 

14 [3, 25] 

86 

Fuhrmann et al. (2014) 

3C 345 

31 [ 20 , 60] 

43 

Schinzel et al. (2012) 

1730-130 

29 [3, 55] 

86 

Fuhrmann et al. (2014) 

3C 454.3 

8 [-4, 20] 

86 

Fuhrmann et al. (2014) 


0 

230 

Wehrle et al. (2012) 


7.2 Start times of radio and 7 -ray activity 

The peaks in the DCF is due to the peak-to-peak difference in the 
radio and 7 -ray light curves. Because of larger emission regions 
and multiple superposed events in the radio light curves, compar- 
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0059+581 
(weekly binned) 



Lag (days) 



Lag (days) 


0059+581 
(monthly binned) 


" S r 

* n t* tt *** 

+ 

ft 

t 

^-1-'-1-'-1-'-1-'-1-'-K 

^-1-'-1-'-1-'-^ 



2009 2009.5 2010 2010.5 2011 2011.5 2012 2012.5 2013 2013.5 

Epoch (years) 



Lag (days) 



Lag (days) 


Figure 6. Light curves for 0059+581 at radio along with weekly (left column) and monthly binned (right column) 7 rays are shown 
in the top panel. Bayesian block representation is plotted in red line oyer the light curve. It is scaled relative to the light curve. 
The DCF for weekly and monthly binned 7 -ray data with significance test using light curve simulation and mixed source method ai'e 
shown in the middle and bottom panel. The Icr, 2cr, and 3cr significance levels are plotted in red, green, and blue dotted lines. 
(The complete figure set is available from the online journal.) 
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ing the delays between the peaks at both wavebands alone might 
be insufficient. Due to this issue and that owing to light-travel de¬ 
lay (Section 7.4), comparing the start times of the activity in ra¬ 
dio and 7 rays can help us constrain the high-energy emission site 
and tell us more about the emission mechanisms at both wavebands 
(Lahteenmaki & Valtaoja 2003; Leon-Tavares et al. 2011). 

Hence, to estimate the beginning of the activity and charac¬ 
terise the variability in the sources at both wavebands, we imple¬ 
mented the Bayesian blocks algorithm which partitions the data 
into piecewise constant blocks by optimising a fitness function 
(Scargle et al. 2013). Using a false-positive threshold of 1% and 
an iterative determination for number of blocks we obtained the 
Bayesian blocks representation for the light curves as shown in 
Fig. 6 . For most sources in 7 rays, we estimated the Bayesian 
blocks using the arrival time of photons extracted from a ROI of 
1° radius centred on the source’s coordinates. Due to the com¬ 
plex structure and variability of the 7 -ray flares in the sources, 
0716-1-714, 1222-1-216, 3C 273, 1510-089 and 3C 454.3, comput¬ 
ing the Bayesian blocks from the photon arrival times was com¬ 
putationally intensive. Hence, we computed their Bayesian blocks 
from the respective light curves. The Bayesian blocks are shown as 
red line in Fig. 6 which in the case of 7 rays are scaled relative to 
the light curve shown as black circles. 

With the aid of the Bayesian blocks, we can notice in most 
of the light curves (see Fig. 6 ) that the start time of an outburst in 
7 rays and in radio to be quasi-simultaneous. From the physical 
scenario, this might imply the emission at both wavebands to be 
associated with the same shocked feature and the seed photons re¬ 
sponsible for the 7 -ray emission might be arising either closer to or 
within the radio core. The former can be constrained to a region up¬ 
stream of the radio core called as the acceleration and collimation 
zone, where the emission feature propagates along a spiral path of 
toroidal magnetic field peaking in the 7 rays as it exits the zone 
(Marscher et al. 2008). Due to synchrotron self-absorption effects, 
this region is opaque at radio frequencies, in which case radio flux 
is in a quiescent state until the moving shock interacts with the 
radio core. The core in the mm-wavebands has the characteristics 
of a standing conical shock that compresses the flow and accel¬ 
erates the electrons. According to the Turbulent Extreme Multi- 
Zone (TEMZ) model of Marscher (2014), a standing shock oriented 
transverse to the jet axis at the vertex of the conical shock can create 
a variable nonthermal seed photon field that is highly blueshifted in 
the frame of the faster jet plasma, leading to rapidly variable 7 -ray 
emission. 

If the seed photons responsible for the 7 -ray flare were to 
arise within the radio core, then the viable explanation would be 
the interaction of a moving shock with the core. This is shown by 
the correspondence of a 7 -ray flare with the ejection of a super¬ 
luminal component from the radio core by various VLBI analyses 
(Jorstad et al. 2001; Schinzel et al. 2012; Jorstad et al. 2013). The 
interaction of the moving shock with the quasi-stationary feature 
downstream of the radio core also contributes to the observed ac¬ 
tivity in the 7 rays (Agudo et al. 2011a). 


7.3 Size of the emission region 

In majority of the correlations, we found the peak of the 7 - 
ray emission in FSRQs to precede those at radio in time-scales 
of days-months. Therefore, to put the correlation results into a 
more physical context, we estimate the distance travelled by the 
emission region, d-,,radio, from the time lags using the relation 


(Pushkarev et al. 2010), 


(5) 


where /lapp is the apparent jet speed in units of speed of light c, 0 
is the jet viewing angle, is the time lag in the observer’s 

frame and 2 is the redshift. We obtained /3app from VLBA moni¬ 
toring of AGNs at 15 GHz (Lister et al. 2013) and the jet viewing 
angle was estimated using the /3app and a variability Doppler fac¬ 
tor as shown in Hovatta et al. (2009). Eor 0235-1-164 we used the 
value of /3app from Agudo et al. (2011b). We were unable to esti¬ 
mate d-).,radio for six sources in Table 2 owing to the lack of /3app or 
variability Doppler factor. The values of d.y,radio for sources having 
/3app and 0 are given in Table 2, along with its lower and upper lim¬ 
its in square brackets. We also estimated d.y,radio for the time lags 
obtained from stacking analysis, by taking /3app and d as a mean 
from the corresponding sample type (see Table 3). 

The distance, d.y_radio, of 7 pc estimated from the stacking 
analysis, corresponds to a projected distance of ~ 0.7 pc for an 
averaged redshift of 0.9 and viewing angle of 3°.3 for sources 
with significant correlation. This corresponds to a projected size 
of ~0.08 mas. Jorstad et al. (2001) estimated the sizes of 43 GHz 
VLBA core for various blazars using multi-epoch observations. 
They have 47 observations for the sources in our sample show¬ 
ing significant correlation, with the average size of 0.1 ± 0.02 mas. 
This when compared with our estimate of 0.08 mas allows us to 
constrain the 7 -ray emission site within the radio core, which is in 
line with the far-dissipation scenario in most cases. 


7.4 Light-travel argument 

The variability in the emission of blazar jets based on the inference 
discussed above, is due to the shock-shock interaction over a fi¬ 
nite size or a time interval. In either case, the radiative cooling time 
of electrons producing the emission in radio through synchrotron 
mechanism is relatively longer than the IC mechanism producing 
the 7 rays. Hence, the decay time-scales are significantly longer 
in radio than in 7 rays. Due to the size/duration of the shock in¬ 
teractions or the evolution of an internal shock, observations are 
affected by the light-travel delay so that variations faster than the 
light-travel time will be spatially unresolved (e.g., Sokolov et al. 
2004; Chen et al. 2011). 

Using light-travel delay argument, Nalewajko et al. (2014) 
pointed out that the temporal coincidence of radio and 7 -ray flares 
alone cannot be used to constrain the site of the 7 -ray emission, if 
the delay between the 7 -ray and radio emission is long enough. 
Based on the assumption that a 7 -ray flare is observed 
when a moving shock interacts with the standing shock, the au¬ 
thors pinpoint the time when the 7 -ray photons were emitted as 
tj,em = t^,obs — (fcore — rj)/c, where Tcore and Tj are the dis¬ 
tance to the radio core and the 7 -ray emission site from the central 
engine and c is the speed of light. Thus, for the 7 -ray flare to be 
produced within the radio core, Nalewajko et al. (2014) proposes 
the relation, (rcore — r-y) <C 2r^ni/3mmCfmm to be satisfied, where 
/SmmCfmm cotTesponds to the size of the component. Here mm cor¬ 
responds to millimetre wavelength. These relations, however, hold 
only when a long-enough time delay is observed. We note that in 
many cases the 7 -ray and radio activity begin simultaneously (see 
Section 7.2), in which case the relation does not hold. 

It should also be noted that the time-scales of variability reflect 
the size of the emitting region and not its distance from the central 
engine. The short time-scales of high-energy variability observed 
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Figure 7. Scatter and boxplots for radio vs 7 -rays for correlated and uncorrelated sources with number of Bayesian blocks (left) and averaged fluxes for the 
5 yr duration (right). The box comprises 75% of the distribution of data with plus symbols denoting the outliers in the boxplots. The median of the distribution 
is denoted by a solid line in the box. The scaling of the boxplots are similar to the scatter plot for both left-hand and right-hand panels. 


in many sources can still be reconciled with the emission region lo¬ 
cated parsecs away from the black hole due to very compact emis¬ 
sion regions embedded within the jet (e.g., Ghisellini & Tavecchio 
2008; Giannios et al. 2009) or due to turbulence in the jet flow 
(Marscher 2014). An alternative possibility is the formation of 
a small emitting nozzle in the wake of a strong recollimation 
(Marscher et al. 2008; Nalewajko & Sikora 2009). 


7.5 Uncorrelated flares 

In more than 50% of the sources, we find no correlation between 
radio and 7 rays. To study the possible differences between the 
correlated and uncorrelated sources, we compared the number of 
Bayesian blocks (indication of flaring activity, left-hand panel of 
Fig. 7) in the two sets of sources. Based on our false-positive thresh¬ 
old, we can expect at least one source to have more than one block 
by chance. 

In the scatter plot, uncorrelated sources with two or less blocks 
in radio are hampered by the sampling. In the case of three blocks 
for uncorrelated sources in radio, we notice a major 7 -ray flare with 
no radio counterpart or radio flare with no counterpart in 7 rays. 
Despite the existing uncertainty, these so-called “orphan” flares 
have been attributed to hadronic processes in the case of former 
(e.g., Bbttcher 2007), while for the latter a complex shock structure 
and possible shock interactions are required for producing the high- 
energy activity (e.g., Aller et al. 2014). According to the TEMZ 
model (Marscher 2014), orphan flares are due to the combined ef¬ 
fects of the modulation in the magnetic field and electron energy 
distribution across different turbulent cells and light-travel delay. In 
uncorrelated sources where four or more blocks are shown in the ra¬ 
dio band, the lack of correlation is most likely due to the rapid vari¬ 
ability in the 7 -ray band. This, especially, applies to the sources, 
0716-1-714, 1510-089 and 3C 279, comprising [10,10], [6,9] and 
[4,18] ([radio, 7 -ray]) Bayesian blocks, which in turn, might imply 
the existence of different IC mechanisms for the 7 -ray flares and/or 
to the presence of multiple emission features which manifests as 
smooth increase in the radio light curves due to the superposition 
of various events, while exhibiting multiple flares in 7 rays. As an 
example for an uncorrelated source, the light curve of 3C 279 along 


with its DCF is shown in Fig. 8 . We performed a two-dimensional 
two-sample KS-test (Peacock 1983), to the estimated number of 
Bayesian blocks in radio and 7 rays for the correlated with the un¬ 
correlated sources. The probability that the distributions of corre¬ 
lated and uncorrelated sources come from the same population is 
less than 4%. 

We also compared the averaged flux for the 5 yr period at 
both wavebands for both correlated and uncorrelated sources. This 
is shown using the scatter and boxplots in logarithmic scale in 
the right-hand panel of Fig. 7. The averaged fluxes show that 
the correlation is significant for brighter sources at both wave¬ 
bands. We also quantify this statement based on a two-dimensional 
KS-test as performed above, which shows that the distribution of 
fluxes in radio and 7 rays for correlated sources to be signifi¬ 
cantly different from the uncorrelated ones with chance probabil¬ 
ity of less than 3%. These results indicate that only the strongest 
flares in the two bands are correlated, as was already suggested by 
Lahteenmaki & Valtaoja (2003). It is possible that we obtain statis¬ 
tically significant results only for the strongest flares, and therefore 
we cannot draw strong conclusions based on the weak sources. 


8 SUMMARY AND CONCLUSION 

The correlation analysis using 5 yr of FennifLAT and 37 GHz ra¬ 
dio data for 55 LAT-detected blazars are presented in this work. The 
cross-correlation analysis of individual sources revealed 16 sources 
from weekly binned and 21 sources from monthly binned to be sig¬ 
nificant at 2.36 ct level. For 15 sources this is the first time a 
significant correlation is reported. In majority of the correlations, 
we find the peaks of the 7 -ray emission to precede those at radio. 
We also stacked the correlations for the whole sample and for 40 
FSRQs and 14 BL Lacs to obtain a significant result and an aver¬ 
age estimate for the corresponding sample type. The time lags for 
the weekly and monthly binned for the whole sample are 80 and 
120 d, corresponding to 47 and 70 d in the source frame. There is 
no significant difference in the time lags of FSRQs and BL Lacs 
upon comparison of weekly and monthly binned results. 

The distance travelled by the emission region was calculated 
from the time lags obtained for the significant correlations. For the 
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Figure 8 . Left: Weekly binned 7 -ray (top) and radio (bottom) light curve of 3C 279. Right: DCF of the light curve with no signihcant correlation. 


whole sample we obtained the distance between the occurrence 
of 7 -ray flare and the peak of the radio flare to be ~ 7 pc (de- 
projected). Two sources - 1633-1-382 and 3C 345 - showed a posi¬ 
tive time lag with radio leading the 7 rays. This suggests that the 7 - 
ray emission could come from downstream of the radio core (e.g., 
Leon-Tavares et al. 2011). 

Using the distance travelled by the emission region (projected 
distance ~ 0.7 pc) we obtained a size of ~ 0.08 mas. We com¬ 
pared this estimate with the average size of the radio core of 
0.1 ± 0.02 mas obtained from VLBA observations as presented 
in Jorstad et al. (2001), thereby allowing us to constrain the 7 -ray 
emission site to be co-spatial with the radio core. 

Bayesian block analysis of the light curves at both wavebands 
shows that sources with significant correlations are more variable 
than the uncorrelated ones. In most of the cases, the blocks also 
shows that the start time of a 7 -ray event corresponds closely to 
those in the radio, implying a co-spatial origin for the 7 -ray and 
radio emission regions (e.g., Agudo et al. 2011a). 

Sources with no significant correlation were compared by us¬ 
ing the number of Bayesian blocks and the averaged fluxes. We 
find that sources with two or fewer blocks were affected by the 
sampling of the light curves, and for sources with four or more 
blocks, the uncorrelation to be due to rapid variability in the 7 
rays. The average fluxes for the correlated sources were higher 
than for the uncorrelated ones, implying that only the strongest 
7 -ray flares are correlated with the radio events, as suggested by 
Lahteenmaki & Valtaoja (2003). 

These results are in favour of the far-dissipation scenario, sug¬ 
gesting that the origin of the seed photons for the high-energy emis¬ 
sion is within the jet. The co-spatiality of the high-energy emission 
region with the radio core is the preferred scenario given the size 
of the emission region. This could imply that they are generated by 
the SSC mechanism, or by the EC mechanism with the seed pho¬ 
tons coming from a sheath layer in the jet or an outflowing BLR. 
Modelling the SEDs using simultaneous observations are needed to 
probe this further. 

Our correlation results are in good agreement with the re¬ 
cent findings of Fuhrmann et al. (2014). Although the effect of syn¬ 
chrotron self-absorbed opacity is unavoidable at radio frequencies, 
we could still place constrains on the location of the high-energy 
emission. 
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